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We investigate how the accuracy and stability of numerical relativity simulations of ID colliding 
plane waves depends on choices of equation formulations, gauge conditions, boundary conditions, 
and numerical methods, all in the context of a first-order 3+1 approach to the Einstein equations, 
with basic variables some combination of first derivatives of the spatial metric and components of the 
extrinsic curvature tensor. Hyperbolic schemes, specifically variations on schemes proposed by Bona 
and Masso and Anderson and York, are compared with variations of the Arnowitt-Deser-Misner for- 
mulation. Modifications of the three basic schemes include raising one index in the metric derivative 
and extrinsic curvature variables and adding a multiple of the energy constraint to the extrinsic 
curvature evolution equations. Redundant variables in the Bona-Masso formulation may be reset 
frequently or allowed to evolve freely. Gauge conditions which simplify the dynamical structure 
of the system are imposed during each time step, but the lapse and shift are reset periodically to 
control the evolution of the spacetime slicing and the longitudinal part of the metric. We show that 
physically correct boundary conditions, satisfying the energy and momentum constraint equations, 
generically require the presence of some ingoing eigenmodes of the characteristic matrix. Numer- 
ical methods are developed for the hyperbolic systems based on decomposing flux differences into 
linear combinations of eigenvectors of the characteristic matrix. These methods are shown to be 
second-order accurate, and in practice second-order convergent, for smooth solutions, even when the 
eigenvectors and eigenvalues of the characteristic matrix are spatially varying. 



I. INTRODUCTION 



The goal of projects such as the ground-based Laser Interferometric Gravitational Wave Observatory (LIGO) and 
the space-based Laser Interferometer Space Antenna (LISA) is to detect gravitational waves, and to use them as a new 
observational window for relativistic astrophysics. A primary source for these gravitational waves is the coalescence 
of binary black holes JO . The highly nonlinear and dynamical merger phase of this coalescence process can only be 
calculated by numerical relativity, and obtaining merger gravitational waveforms, both for theoretical understanding 
and for detection, is dependent on long-term stable and accurate numerical evolutions. A worldwide collaboration 
of numerical relativists, physicists, mathematicians, and computer scientists has devoted considerable effort over the 
last 20 years to develop 3D codes to calculate black hole merger gravitational waveforms, and significant progress 
has been made, especially in the last few years. However, more groundwork is required before calculations of 3D 
binary black hole merger templates for a variety of scenarios can be completed. Greater understanding of equation 
formulations, boundary conditions, and dynamic gauge conditions, and the use of advanced numerical methods, is 
essential to achieve this goal. We believe that an important foundation for this understanding is extensive testing and 
analysis in ID and 2D. Choptuik's discovery of black hole critical phenomena in spherically symmetric gravitational 
collapse is an example of the potential of careful numerical work in ID. 

This paper reports the methodology, results, and analysis of calculations of ID nonlinear colliding gravitational 
planewave spacetimes. We have chosen to investigate hyperbolic formulations of the Einstein equations, as they are 
well-posed, they can be treated with advanced numerical methods, and they can help in the analysis of boundary 
conditions ||,|| . We call a set of equations hyperbolic if the characteristic matrix can be diagonalized with a complete 
set of eigenvectors and real eigenvalues, following LeVeque This is called strongly hyperbolic || in much of 
the literature. The lapse and the shift are evolved during each time step in a manner which is consistent with a 
simple hyperbolic scheme. Between time steps, the lapse and the shift are reset according to conditions which are 
unconstrained by the need to preserve hyperbolicity. In this way, the evolution of the hypersurfaces and spatial 
coordinates can be controlled to prevent large gradients, coordinate pathologies, and instabilities. Some of the 
redundant variables of a hyperbolic formulation can also be reset between time steps. This resetting can have positive 
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or negative effects on accuracy and stability, depending on the eigenmode structure of the reset system. Finally, we 
find that boundary conditions should not be based naively on the eigenmodes of the hyperbolic decomposition for 
two reasons: (a) satisfying the constraint equations at the boundaries generically requires the presence of incoming 
eigenmodes, and (b) even whether the "physical" eigenmodes are purely outgoing at the boundaries is gauge dependent. 

Many ways of formulating evolution equations for the spatial metric in Einstein's theory of General Relativity are 
possible. The most thoroughly tested formulation in numerical relativity is the Arnowitt-Deser-Misner (ADM) set of 



equations Q| . The standard ADM equations in vacuum are 

(a t -Cp)h ij =-2aK ij , (1) 

{d t ~ C p )K l3 = -a {ij + a[^R l3 + KK l3 - 2K i l K lj ], (2) 

(d t - C p )K^ = -cP j + a[ {3) R l j + KK'j}. (3) 



In these equations, 2£y is the extrinsic curvature, K = K l i, j3 l is the shift, a is the lapse, h%j is the 3-metric, and ^'Rij 
is the 3D Ricci tensor. The vertical bar represents a covariant derivative taken with respect to the 3-geometry. Eq. 

evolves what we call the "mixed" form of the extrinsic curvature tensor. The energy and momentum constraint 
equations are, respectively, 

£ = 1/2[<®R - KjlTj + K 2 } = 0, (4) 

M i = K i j lj -K li = 0. (5) 

While successful calculations using the ADM formulation have been done in 2D, 3D calculations generally crash after 
just a few dynamical times. 

Alternative formalisms include many versions of hyperbolic systems, which add redundant variables and/or con- 
straint terms to the equations to allow a complete set of eigenmodes describing evolution along the characteristics. As 
indicated by Reula there are an infinite number of hyperbolic formulations. We focus on variations of relatively 
simple schemes proposed by Bona-Masso (BM) (8| and Anderson- York (AY) || , in which the characteristics propagate 
either at local light speed or along the hypersurface normals, and in which the variables include first derivatives of 
the metric. 

Initial attempts at using hyperbolic methods in 3D were based on the BM formulation p| , but did not use numerical 
methods which take advantage of the eigenfields of the system. These codes were not much more successful than ADM. 
Non-hyperbolic Baumgarte-Shapiro-Shibata-Nakamura (BSSN) schemes M,0, based on conformal decomposition 
of the metric, have shown considerable success in improving the stability of 3D calculations for weak and strong 
gravitational fields and a variety of spacetime slicings fl2[| . Alcubierre et al. jl3| report that a BSSN scheme, combined 
with excision and certain dynamic gauge conditions, allows accurate numerical evolutions of 3D distorted dynamic 
black holes up to hundreds of dynamical times. 

In the context of considering only first derivative variables, a great variety of hyperbolic schemes have been proposed 
that involve adding constraint terms to the equations |]l4]^19||. Kidder, Scheel, and Teukolsky |nj examine a rather 
general class of such schemes, which include the AY || and Frittelli-Reula |Q formulations as special cases. Among 
these schemes are some which allow for long-term evolution of a Schwarschild black hole in 3D. 

In this paper, we explore ways of using hyperbolic methods that combine superior accuracy with gauge conditions 
which maintain stability at least for the limited dynamical times we can explore with plane waves. Three basic 
first order systems are studied: BM, AY, and ADM. Hyperbolicity is obtained in our BM and AY formulations by 
adding momentum constraint terms to the ADM equations, as in the standard formulations. The BM, AY, and 
ADM formulations are modified by using "mixed" forms (with one index raised) of the extrinsic curvature and metric 
derivatives as variables. The BM formulations are further modified by resetting redundant variables, which gives an 
overall ADM-like evolution. Further, all our formulations are varied by adding a multiple of the energy constraint 
to the evolution equations for the extrinsic curvature. Specifically, —na£hij/2 is added to Eq. (Q), and —na£S l j/2 
is added to Eq. (||), where n, the energy constraint coefficient, is an arbitrary real number. The ADM formulation 
is actually hyperbolic as long as the longitudinal-transverse components of the metric and extrinsic curvature can be 
assumed to vanish identically, and n<0or0<n< 1. Comparisons of results from our various ADM, BM, and 
AY calculations allow us to identify and analyze aspects of equation formulation which significantly improve accuracy 
and/or stability. These are mixed variables, a separation of the constraint error speeds from the other characteristic 
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speeds of the system, and maintaining long-term effective hyperbolicity (taking into account resetting of redundant 
variables, but ignoring deviation from strict hyperbolicity due to resetting the lapse and the shift). 

While the same energy constraint terms as specified above are present in the standard BM formulation, numerical 
implementations have been carried out, as far as we are aware, only for n = (the Ricci evolution system) and 
n = 1 (the Einstein evolution system). Adding these energy constraint terms to the AY formulation is a special case 
of the more general Kidder, Scheel, and Teukolsky Jl9| schemes. Shinkai and Yoneda Jl5|-|l7j analyzed the stability 
and accuracy properties of first order hyperbolic systems using Ashtekar's connection variables in plane-symmetric 
spacetimes, and found that the addition of multiples of the constraints to the dynamical equations improved accuracy 
and stability. These results have been extended to ADM systems of equations Jig ). However, the approach of Shinkai 
and Yoneda is different from ours in that they consider the constraints as independent dynamical variables. 

Gauge choices in most previous implementations of hyperbolic formulations have been limited in order to preserve 
the hyperbolicity of the system. Since no time derivatives of the lapse and the shift occur in the dynamical equations 
for the other variables, the lapse and the shift can be reset arbitrarily at any time during the numerical evolution, 
as pointed out by Balakrishna et al. pQ] . Our gauge evolution maintains strict hyperbolicity during each time step, 
but the lapse and shift are reset periodically between time steps in order to control the long-term evolution of the 
coordinate system. The resetting may be accomplished by imposing algebraic conditions, by solving elliptic equations, 
or by evolving the lapse and/or shift through dynamical equations implemented independently of the main hyperbolic 
system. 

Poor boundary conditions can result in the introduction of instabilities or inaccuracies into the numerical grid. 
In numerical relativity, boundary conditions have usually been rather crudely implemented. Some sort of outgoing 
radiation conditions are imposed on all components of the metric, or boundary conditions are based on an analytic 
exterior solution |l^,|2l| . One attraction of hyperbolic methods has been the possibility of basing boundary conditions 
on the eigenmodes of the characteristic matrix. However, it is clear from our planewave calculations that, particularly 
for the "non-physical" eigenmodes involving the non-transverse-traceless parts of the metric, making the amplitudes 
of the incoming eigenmodes at the boundaries zero can lead to serious violations of the energy and momentum 
constraints. Furthermore, what constitutes an incoming eigenmode is dependent on the formulation of the equations 
as well as on gauge conditions. Even imposing purely outgoing boundary conditions on the "physical" eigenmodes 
of the hyperbolic system is not strictly correct, as nonlinear coupling between the "physical" and "non-physical" 
eigenmodes in the source terms can generate a gauge-dependent admixture of outgoing and incoming "physical" 
eigenmodes. Our boundary conditions are based on quadratic extrapolation of the variables from inside the grid to 
the first ghost cells on either side of the grid. The ghost cell values are then corrected to make sure the constraint 
equations are satisfied on the boundaries. For ID plane waves, projection of the Weyl tensor onto a null tetrad gives 
a gauge-independent measure of the left and right-going components of the gravitational radiation. Our numerical 
solutions for colliding plane waves show that as the wave packets leave the grid, the incoming components of the Weyl 
tensor are in fact zero even though there are non-zero incoming "physical" eigenmodes of the characteristic matrix. 

Our focus in applying hyperbolic methods to the Einstein equations is on achieving second order accuracy for smooth 
solutions, when the eigenvectors and eigenvalues of the system are a function of position. Finite difference methods 
such as MacCormack, Lax-Wendroff, and staggered leapfrog p^ |, which are often used in numerical relativity, give 
good second order accuracy for smooth solutions, but standard wave propagation algorithms for hyperbolic systems as 
presented by LeVeque p3| are not second order accurate when the eigenvectors and eigenvalues are spatially varying. 
LeVeque suggested a new wave propagation method for variable coefficient flux problems which we develop and 
apply to our ID nonlinear gravitational planewave calculations. We show in the Appendix that the new methods are 
formally second order accurate even with varying eigenvectors and eigenvalues, and verify second order convergence 
in our numerical results. 



II. EVOLUTION EQUATIONS 

The most general spatial metric for a nonlinear ID plane wave traveling in the re-direction is 

ds 2 = h xx dx 2 + hyydy 2 + h zz dz 2 + 2h yz dydz, (6) 

in which h XXl h vv , h zz , and h yz are functions of x alone. We will restrict our discussion to a diagonal metric in this 
paper. 

The standard ADM evolution equations are first order in time and second order in space. Most hyperbolic formalisms 
are first order in space and time, and incorporate first derivatives of the spatial metric as additional variables. The 
derivative variables as defined by BM are 
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In our applications with a diagonal metric, we find that switching to the "mixed" variables 

D k j = \{d k h a )h l \ Kf =K a h l i 



(7) 



(8) 



from the "lowered" variables, D^ij and Kij, improves accuracy significantly without other complications. However, 
with a non-diagonal metric, Dki 3 and are not symmetric in i and j, and the evolution equations for Dm 3 acquire 
complicated source terms. 

Below, we present the first order evolution equations for ID plane waves travelling along the x-direction and 
described by a diagonal spatial metric, using our mixed variables. The equations with D k ij and Kij as variables are 
given in the BM papers Q. A few points need to be made about notation. First, since our ID problem involves 
derivatives only in the x-direction, we simplify our notation D x i 3 — > Di J . Second, a prime indicates a spatial 
derivative with respect to x. Third, our symbol for the shift is simply /3 instead of (3 l . We suppress the index on the 
shift because there is only one non-zero component in this ID case. 

The evolution equations for hij are obtained from the definition of the extrinsic curvature of the hypersurfaces, Eq. 
(Q), and are 



d t h xx = 2h xx [f3D x x + p' - aK x ° 



(9) 



dthyy = 2h yy [pD y v - aK y y ], 8 t h zz = 2h zz \pD z z - aK z 



(10) 



The evolution equations for ZV are obtained by taking the time derivative of D^i 3 in Eq. (H) , and interchanging 



space and time derivatives. The resulting equations are 

d t D x x + d x [-PD x x -f3' + aK x x 



0. 



(11) 



d t D y y + d x [-[3D y y + aK y y ] = 0, d t D z z + d x [-(3D z z + aK z z ] = 0. 



(12) 



The Ki 3 variables are evolved from the Einstein equations, Eq. (|3|). We include the addition of an arbitrary multiple, 
n, of the energy constraint in these equations. After organization into a conservation law form, the evolution 
equations are 



3tK x * + a 

-f3'K x x + a 



a a 



-/3K X X + —(- + D y y + D z 

xx V a 

1 



K x x Ki 



hxx V a 
1 /'a' 

w 



y y z 2 



D x 



D y y + D z 



-o£, (13) 



d t K y v + d x 



-p Ky y + —D y y 



-0K v y 



a 



(14a) 



d t K z z + d x 



-S3K Z Z + ^D z 



= -0K Z Z + a 



n 



where we write a£ so that the division between the flux terms and the source terms is apparent: 
o£ = —d x 



( Dy y + D z z ) 



a[K x x (K y y + K z z ) + K y yK z z } 



~-Dx X ) (D y y + D z z ) - (D y v D y y + D y yD z z + D Z Z D Z 



(14b) 



(15) 
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III. GAUGE EVOLUTION 



We let the lapse and the shift evolve during each time step according to a prescription which simplifies the hyperbolic 
system, and we periodically reset the lapse and the shift between time steps to control the longer term evolution of 
the coordinates and to keep gauge pathologies from developing. We defer discussion of resetting gauge conditions to 
Sec. 0. Here, we discuss how the gauge evolves between resettings. 

For a hyperbolic formulation of the equations, the natural choice for the lapse between resettings is the Choquet- 
Bruhat algebraic gauge condition p5|pr| |, because it simplifies the fluxes and source terms in the hyperbolic system 
of equations considerably. This gauge condition is 

a = Q^det{h l3 ), (16) 

where Q is a specified function of x, t. 

We vary the Choquet-Bruhat algebraic gauge condition by making Q and Q' (which equals d x Q) variables in the 
hyperbolic system, rather than specified functions of x, t. We choose Q and Q' as variables so that Q' can be included 
in the flux of K x x as part of the hyperbolic system. Otherwise, Q" would have to be considered part of the source of 
K x x , and evaluating Q" from the lapse involves second derivatives of the metric. By advecting Q and Q' along the 
hypersurface normals, we incorporate them into the hyperbolic system in a consistent way. Our advection equations 
are 

d t Q-pQ' = Q, d t D Q -d x [/3D Q }=0, (17) 

where Dq = Q' /Q. Our advection of Q corresponds to harmonic slicing p5[ . 

There is a danger with resetting the lapse and the shift, in that fluctuations in 0' and Dq can feed back on one 
another through the evolution equations for D x x and K x x . The resetting gauge conditions of Sec. |v| imply that a 
fluctuation in K x x is balanced by a fluctuation in /3' , and a fluctuation in D x x is balanced by a fluctuation in Dq. 
For certain time intervals between resetting, if these fluctuations propagate at different speeds, they may drift in such 
a way that they reinforce rather than cancel over much of the time interval. Although the standard procedure is to 
keep the shift constant in hyperbolic schemes, we find that if we advect Dq with f3' constant, such a positive feedback 
can occur, resulting in a runaway instability. However, if we advect both Dq and [3' along hypersurface normals, the 
evolution is stable. Our advection equations for and [3' are 

9t0-pp = O, dtp - d x [(30'} = 0. (18) 



IV. CONSTRAINT EQUATIONS 

The energy and momentum constraints must be satisfied by the initial conditions and throughout the evolution. 
We use these constraints to obtain the initial conditions. We do not impose the constraints during the evolution 
of the dynamical equations. However, we do insure that the boundary conditions are consistent with the constraint 
equations, and we use the constraints to check for accuracy and convergence as the numerical evolution proceeds. The 
energy and momentum constraint equations are, respectively, 



£ 



-d x 



1 



W + D z z ) 



1 



[Dy V D y v + D y y D z z + D / D / + D X '(D V V + D / )} + K x x (K y y + K/) + K y y K z 



(19) 



M x 



-d x (K y y + K z z ) 



D V K v 

Uy JXy 



D Z Z K Z 



D Z Z )K X = 0. 



(20) 



V. RESETTING GAUGE CONDITIONS 



The lapse and shift are periodically reset between time steps in order to implement a dynamic spacetime slicing 
which is unconstrained by the need to maintain a hyperbolic system. Our resetting gauge conditions are chosen to 
prevent pathologies and/or strong gradients from developing in the hypersurfaces and spatial coordinates, and to help 
stability properties at the boundaries of the grid. 
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Changes in spacetime slicing which maintain the explicit planar symmetry only directly impact K x x . Nonlinear 
source terms in the evolution equation for K x x have the potential to generate runaway growth of K x x when K x x is 
positive. Our lapse resetting condition drives K x x toward a small negative value to insure against this. In addition, 
a negative K x x implies that the proper distance between hypersurface normals displaced in the x-direction increases 
with time. Together with our shift resetting condition, which keeps h xx roughly constant, this results in hypersurface 
normals which point outward at the boundaries of the computational domain. Some features in K x x and D x x 
potentially associated with instability advect along hypersurface normals (see Sec. IX), and are then advected out of 
the grid before they can do much harm. 

The equation for the lapse is derived by imposing the condition at the time of resetting 

d t K x x - (3d x K x x = -T[K X X - {K X X ) T ], (21) 

where (K x x )t is a specified "target" value, and T is a damping constant, which is chosen to be comparable to the 
characteristic frequency of the waves we are propagating. Substituting this condition into the evolution equation for 
K x x (Eq. ([l3|)) and simplifying using the energy constraint, we obtain our lapse resetting condition, 



a. 



a' 



T(K X X - {K X X ) T ) + K X X K X - K*K* + 



/>/(£-)■ <->->> 



To limit initial transients in the lapse, given our initial condition K x x = 0, the target value is made proportional to 
(1 - e- r */ 4 ). 

In our colliding wave calculations, Eq. (|2|) as it stands can cause the lapse to become negative at the edges of 
the grid, if the second derivative of the lapse becomes too negative. To prevent this, we replace S, the expression in 
square brackets in Eq. (p2|), by 

S — S (23) 

y/l + (S/S Hm ) 2 

when S is negative, so S > — \Su m \. A side effect of the limiter is to allow K x x to become more negative than its 
target value. 

We choose an equation for the shift so that at the time of resetting, h xx is advected along hypersurface normals: 

d t h xx - (3d x h xx = 0. (24) 
Substituting this requirement into the evolution equation for h xx (Eq. @), we obtain the shift resetting condition, 

8 x p = aK x x . (25) 



VI. HYPERBOLIC SYSTEMS 

The evolution equations presented in Sec. |n| have been cast in a first order, flux-conservative form, represented by 
the following set of I equations 

d t q + cUF(q)] =S(q), (26) 

where q is a vector of / variables. The flux vector is given by 

F(q) = A(x)q, (27) 

where the / x I characteristic matrix A(x) is the flux Jacobian, 9 q [F(q)]. The system is hyperbolic if A(x) has a 
complete set of eigenvectors and real eigenvalues. 
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A. Modified Bona-Masso Formulation 



The standard BM formulation || creates a hyperbolic scheme by introducing the redundant variables Vi, which are 
defined as 



Vi — Dik k — D k ki 



V x = D y v + D z 



The momentum constraint is used to evolve V x : 

d t V x + d x [-l3V x ] = a[D y y K y y + D Z Z K Z Z - (D y v + D Z Z )K X X - (a'/a)(K v y + K z 



(28) 



(29) 



We densitize the lapse according to the Choquet-Bruhat algebraic condition (see Sec. IB), which simplifies the 



standard BM system of equations considerably. With a — > Q \/ det(hij) and (A x = d x In a) — ► Dq+D x x +D y v +D Z Z , 
the fluxes for AV reduce to 



F{K X 



-PK X X 



D Q + D X X + (2-^)v x ) 



(30) 



F(K y y) 



h xx 



n y 



V x ), F(K Z Z ) 



-f3K z 



h x 



-( 



D. 



-V, 



(31) 



Our advection of Q and Dq, as described in Sec. BI, corresponds to harmonic slicing, a special case of the standard 
BM lapse evolution equation. We also advect (3 and j3' , whereas the standard BM formulation specifies the shift as a 
known function of x and t. 



B. Modified Anderson- York Formulation 



The AY formulation differs from the BM scheme in how the momentum constraint is used to make the system 
hyperbolic. The AY scheme eliminates the need for the BM redundant variables Vi by incorporating the momentum 
constraint into the evolution equation for the fkij variables, which are defined as 

fkij = Dkij + hkiVj + hkjVi. (32) 

The Vi variables in this equation are not separate variables, but rather denote the combinations of the _D's given in 

The AY formulation replaces the BM D^j with fkij , which are simply the spatial metric derivative terms in the K^j 
fluxes of the BM formulation. Using this as a guide, we generalize the AY scheme (whose original form is restricted 
to Ricci evolution, n = 0) to allow for non-zero energy constraint contributions. This leads to 

TL 

fkij = Dkij + h^Vj + hkjVi — —Vkhij. (33) 



The generalization in Eq. (33) works as long as the inverse transformation from fkij to Dkij exists, which is the case 
for n 1. An evolution equation is obtained for fkij from Eq. ( [33] ) by using the momentum constraint to eliminate 
the time derivative of the Vi variables. For our modified AY scheme, we then raise one index so that fk% — fkuh 1 ^ 
are our basic variables. A hyperbolic system results without the need for the BM redundant variables Vi. 
For the diagonal metric planewave case under consideration, Eq. (B3J) reduces to 



U* =D x x +(2-^)V x , (34) 



fy V = D y y -^V x , f z z =D z z --V x . (35) 

We have simplified our notation in that fki — > f% for this ID problem. Notice that f yx v = V x and f zx z = V x , which 
contribute to fluxes in the y and z directions, are not zero. Bowever, with planar symmetry the divergence of these 
flux components vanishes identically. 
The evolution equations for are 



7 



d t f x x + d x [~Ph x + aK x x ] - (2 - |) aC, 



(36) 



71 ft 

dtf y y + d x [-(3f y y + aK y y ] = --aC, d t f z z + d x [-[3f z z + aK z z ] = --aC, (37) 



where 



C = [DyVKy* + D/-K/ - (D y y + D Z Z )K X X - (a / /a)(lf/ + K/)]. (38) 
The .D's in Eq. (p8[) are not separate variables, but denote: 



o n 

DS^fS-lj-l) [//+//], (39) 



Z) "" = 2(1^0 ~ + nf/]l D/ = W^n) [{2 " )L ^ " 1 m 



These relations are the inverse transformation of the system of Eqs. ([34|) to One can see that n = 1 is not 

allowed. 

The Ki 3 evolution equations are the same as in our modified BM scheme, with the understanding again that the 
D's in the source terms are not separate variables, but the above linear combinations of /'s (Eqs. (|39|) to (ff(])). The 
fluxes are defined so the fi° variables can replace the expressions involving ZV in the fluxes of our modified BM 
scheme. The following fluxes result: 

F(K x *) = -pK x x +-?-(D Q + f x x ), (41) 

f^XX 

F(K y y) = -f3K y y + -^f y y, f{k/) = -px/ + -£-f x *. (42) 

f^xx flxx 

The AY formalism imposes the Choquet-Bruhat algebraic condition on the lapse, as we did in our modified BM 
scheme. The evolution of the lapse and the shift between gauge resettings is treated in exactly the same way as in 
our modified BM formalism. 



C. Modified Arnowitt-Deser-Misner Formulation 



The simplest of the hyperbolic schemes we present is our modified ADM formulation, which consists of Eqs. (Q) to 
@, (0), and @, with a = Qy^del(h~) and a' /a = D Q + D x x + D y v + D/. This system is hyperbolic when the 
metric is diagonal ifn<0or0<n<l. The fluxes for are 



F(K X X ) = -0K X X + -£L[d q + D x x + (2 - I) W + D z 



(43) 



W)=-w+£-[(i-f) 



2) D * - 2°* 



F(K Z 



-0K Z * + — 



n\ 

2)°' 



2 



'v 



(44) 



The hyperbolicity of our modified ADM system of equations breaks down for n — and n > 1. Although our ADM 
formulation at n = is non- hyperbolic, it is stable. At n — 1, however, the system is both non-hyperbolic and on the 
verge of being unstable. For n > 1, the equations are elliptic, giving unstable exponential growth of errors. 
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D. Wave Modes 



1. BM 



The hyperbolic system of equations obtained from the modified BM formulation described in Sec. VIA is 

a t q + a x [A(a!)q]=S(q), 

where 



(45) 



(D x x 




D y y 




D z z 




K x x 




K y y 




K z z 




v x 




Dq 




\ & 


) 



(46) 



and 



A(x) 



{ 


-p 








a 














-1\ 







-p 








a 






















-p 








Q 













a 








~P 








hx~x~ — i) 


a 
hxx 










a 
hxx 








-P 





a / n \ 
fi„ V 2 / 
















a 
hxx 








-/3 


a / n \ 
hxx V 2 } 




























~P 































-0 





\ 


























-PJ 



(47) 



The nine eigenmodes of the homogeneous system are obtained from the characteristic matrix, A(ir). 
eigenmodes travel along the light cones. They are: 



/hxa 



[D x x + D Q + (2 - f ) V x ] ± 



> speeds 



-P± 



a 



Six of the 



(48) 



The remaining three eigenmodes are simply the variables V x , Dq, and P' , which travel along the hypersurface normals, 
with speeds — p. 

The eigenmodes of the characteristic matrix, however, do not necessarily describe how solutions of the full nonlinear 
system of equations propagate. It is a special property of planewave systems that eigenmodes of the full nonlinear 



purely 



system of equations exist which consist of purely right-going waves with K y v 
left-going waves with K y v ± K z z = —{D y y ± D z z )/\/h xx , and D x x = K x x — 0. These are solutions of the Einstein 
equations in a gauge with a = 1 and P' = 0. In our nonlinear colliding plane wave calculations, our initial conditions 
are such that the waves have this form. The right-going wave is in the left half of the grid, the left-going wave is in 
the right half of the grid, and they are just at the point of colliding. When discussing solutions of the full nonlinear 

) , the constraint 
After the waves 



system of equations, we refer to the transverse-traceless quantities (D y v ~D Z 



quantities (D y v + D z 



± K z z = {D y y ± D z 



and (K y v -K z 



and (K y v + K z z ), and the longitudinal variables D x x j\Jh xx and K x 



pass through each other, it is only approximately true that the transverse-traceless quantities have the form of purely 
right-going and purely left-going waves as described above and it is not at all true that the constraint quantities have 
this form. 

The characteristic speeds apply to small amplitude, short wavelength perturbations in the variables, so that the 
principal terms (which are first derivative terms) dominate over the source terms. The disturbances in the constraint 
quantities which propagate along the characteristics will generally be constraint-violating because the constraints 
explicitly tie the principal terms to the nonlinear source terms, and require that they cancel. The longitudinal variables, 
D x x j \/h xx and K x x , are not eigenmodes of the homogeneous system. In the full nonlinear system, D x x / \Jh xx and K x x 
have some features which propagate along the light cones, and some features which propagate along the hypersurface 
normals. The propagation of the longitudinal variables is strongly dependent on the choice of gauge. 
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2. AY 



There is a complete set of eight eigenmodes of the modified AY homogeneous system of equations. The six 
eigenmodes which travel along the light cones are 



7 j i -(f x * + D Q )±(K x *-£), 

7C ±Ky ' vc ±Kz 1 v 



speeds = -/? ± (49) 



The remaining two eigenmodes are the variables Dq and (3', which travel along the hypersurface normals, with speeds 



3. ADM 



For the modified ADM homogeneous system, the eigenmodes, which form a complete hyperbolic system for n < 
or < n < 1, consist of "longitudinal" and "physical" eigenmodes propagating along the light cone, 



y±= [nD x x + uDq + (2 - f ) (£>„» + D z z )] ± [nK x x - + (2 - f ) (K y y + K/)] , 1 

j-[D y y-D z ']±[K v y-K z '\ J 



"constraint" eigenmodes propagating inside the light cone for < n < 1, 

[£>/ + £>/] ± [K y y + X/] } speeds = -f3± -£= y/T= 



speeds — —(3 ± 



n, 



(50a) 



(50b) 



and the Dq and eigenmodes with speeds — (3. Hyperbolicity fails for n = because the "longitudinal" eigenmodes 
are not independent of the "constraint" eigenmodes, for n = 1 because the two "constraint" eigenmodes are not 
independent of each other, and for n > 1 because the "constraint" eigenvalues are complex. 



VII. BOUNDARY CONDITIONS 



Since in numerical relativity, computations are usually performed on a limited grid within a much larger space, the 
boundary conditions should be designed to be consistent with how waves propagate while they are still inside the 
grid. Even more important, since the evolution equations admit constraint-violating solutions, constraint violations 
will propagate into the grid unless boundary conditions are carefully designed to suppress them. 

Consider the "constraint" eigenmodes. They are [D y y + D z z — nV x ]/\/h xx ± [K y y + K z z ] in BM and AY (though 
expressed in terms of different variables), and y/1 — n [D y v + D z z ]/y/h xx ± [K y y + K z z ] in ADM. Even for the same 
value of the energy constraint coefficient n, what is outgoing in BM and AY is different from what is outgoing in 
ADM. Furthermore, for a given solution, the amplitudes of the BM, AY, and ADM modes depend on n. Whatever the 
correct boundary condition, its effect on the solution should be independent of the equation formulation. The relative 
amount of right and left-going "constraint" modes is also gauge dependent, in the sense that the choice of boundary 
conditions in solving the constraint equations in the initial conditions is a gauge choice, and this affects the relative 
values of (D y y + D z z ) and (K y y + K z z ) at all later times. The initial conditions symmetric about the midpoint of 
the grid at x = 10 give purely incoming "constraint" modes for n = ([D y v + D z z ]/y / h xx — ±[K y y + K z z ] on the 
left /right edges of the grid) initially and at all times until the effects of the wave collision reach the boundaries. 

The "longitudinal" eigenmodes involving D x x and K x x are also formulation dependent, since they are different in 
ADM from what they are in BM and AY, and they depend on n in all three formulations. There is gauge freedom to 
pose any boundary conditions one likes on these modes, but a poor choice might give rise to singularities in D x x or 
K x x inside the grid. 

The "physical" eigenmodes [D y y — D z z ]/ 1 \Jh xx ± [K y y — K z z ] are the same in all three formulations, and are 
independent of n. However, their time evolution is gauge-dependent because the nonlinear source terms in their 
evolution equations involve the gauge-dependent constraint quantities. With our choice of initial gauge, the amplitudes 
of (D y y — D z z )/ \fhx~ x and {K v y — K z z ) differ by about 3 per cent after the wave collision, so there is typically a 3 
per cent admixture of the incoming "physical" cigenmode as the outgoing waves approach the boundaries. 
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A gauge-independent measure of the amplitudes of left and right-going gravitational waves can be obtained by 
projecting the Weyl tensor onto a complex null tetrad, as in the Newman- Penrose spin coefficient formalism [ ^7f , 



*± - R (t)(y)(t)(y) ~ R(t)(z)(t)(z) + R(x)(y)(x)(y) ~ R(x)(z)(x)(z) T 2[R (t)(y)(x](y] - R( t )(z)(x)(z)\ 



(51) 



for right/left propagation. A purely right-going wave would have = 0. Our numerical results indicate that plane 
waves after a collision are indeed purely outgoing by this standard. As an outgoing wave boundary condition, using 
the evolution equations to evaluate the time derivatives of the extrinsic curvature in the Riemann tensor, this becomes 



0,- 



(D y y - d 



1 (£>/ - D z z ) 



l(D y y + D z *) 




T (K y y-K/) 



+ 2 



W~D,') 



T {K y v + K z 
T (K y y - K z 



= 



±[K y y 



(52) 



K z z ] if and only if 



at the right/left boundaries. This expression is consistent with [D y v 
[D y v + D z z }/Vh^ = ±[K y y + K/'}. 

Since conventional outgoing wave boundary conditions are not appropriate, our boundary conditions are based on 
a smooth second order extrapolation of the variables, which is corrected to make sure the energy and momentum 
constraint equations are satisfied on the boundaries. Eq. ( [52| ) could also be imposed at the b oundarie s to further 
improve the extrapolation, but we have not tried this. Our procedure is detailed further in Sec. VIII B 3. 



In addition to the eigenmodes discussed above, there are eigenmodes propagating along the hypersurface normals, 
which can be incoming or outgoing, depending on the sign of the shift on the bound aries . It s eems to be important 
for stability that the hypersurface normals do not point into the grid (see Sees. [X C and IX Ej ). 

Our results show that quadratic extrapolatio n with out correction for the energy and momentum constraints produces 
a significant but not dominant error (see Sec. IX D| ). However, errors from imposing outgoing boundary conditions 
on the "constraint" eigenmodes, or from using standard constant extrapolation, would swamp all other errors as they 
propagate into the grid. Standard constant extrapolation, which gives the same values for the variables, and therefore 
the fluxes, in the ghost cell and adjoining physical cell, also eliminates the incoming "constraint" eigenmodes. 



VIII. NUMERICAL METHODS 



A. Strang Splitting 

As described in Sec. [v| , all of the formulations we tested, both hyperbolic and non- hyperbolic, are in first order, 
flux conservative form. We solve all these sytems of equations using a Strang-split method J22|. In this method, 
the homogeneous transport part of Eq. (|2^) and the contributions from the source terms are treated separately. In 
particular, the following straightforward system of ordinary differential equations is first solved over half a time step 

d t <i = S(q). (53) 

Then, the transport part of Eq. (p6|), which contains the flux terms, is solved over a full time step 

d t q + d x [F(q)} =0. (54) 



Our methods for solving the transport step are discussed in Sec. VIII B below. The calculation is completed by again 
solving Eq. (^3|) over half a time step. 

We choose to use the Strang-split method because it is simpler in the context of how we are handling boundary 
conditions. An iterative scheme such as the MacCormack method [p8| requires repeated implementation of the 
boundary conditions each time step. However, in the Strang-split scheme, the boundary conditions are imposed 
only once during each time step. The fewer applications of the boundary conditions in the Strang-split method 
is advantageous because we are using quadratic extrapolation to obtain ghost cell values. Quadratic extrapolation 
amplifies any jitter at the boundaries, and the frequent application of quadratic extrapolation in iterative schemes 
such as MacCormack could easily lead to an instability. 
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B. Transport Step 



In the transport step, we solve Eq. ( p4[ ) both with a finite difference method and with a wave propagation approach, 
which takes advantage of the eigenfields of a diagonalizable hyperbolic system. Advanced numerical methods for 
diagonalizable hyperbolic systems introduce limiter functions to resolve sharp discontinuities that typically arise in 
hydrodynamics problems. A smooth problem can be solved just as accurately and more efficiently with a finite 
difference method. In vacuum general relativity, discontinuities may or may not arise, depending on the gauge 
conditions. Commonly used gauge conditions lead to steep gradients near black hole horizons. One can deal with these 
gradients by using high resolution methods requiring diagonalizable hyperbolic formulations; or, one can dynamically 
adjust the gauge conditions so as to avoid the steep gradients altogether [ p9| . 

Whether one uses a finite difference method or a sophisticated hyperbolic technique, it is important to have a 
numerical scheme which is fully second order accurate for smooth solutions and generalizable to black hole spacetimes 
and higher dimensions. It is straightforward to devise a finite difference scheme based on a Taylor series expansion 
which is formally second order accurate. High resolution Riemann-based wave propagation algorithms introduced 
by LeVeque |23), which decompose Aq across a grid cell interface into a linear combination of eigenvectors of the 
A(x) matrix, are applicable to a wide variety of diagonalizable hyperbolic problems. Flux differences are calculated 
from the Aq decomposition. We refer to these algorithms as "standard wave decomposition" methods. However, the 
standard wave decomposition methods are not second order accurate for smooth solutions when the characteristic 
matrix A(x) is a function of position, because the changes in A(x) across cell boundaries as well as Aq's must be 
accounted for in flux differences. In numerical relativity problems, A(x) depends on the lapse, the shift, and the 
spatial metric, and can have gradients comparable with the gradients in q. 

LeVeque has suggested a wave propagation approach for solving variable coefficient flux problems based on splitting 
up the jump in F(q) rather than the jump in q [Q. We refer to this approach as "flux-based wave decomposition". 
We develop and apply this method to solve the Einstein equations for ID nonlinear plane waves as described below 



in Sec. VIII B 1 . We show in the Appendix that flux-based wave decomposition methods are formally second order 



accurate for sufficiently smooth solutions for arbitrary smooth variations of the eigenvalues and eigenvectors (see also 
Bale et al. p0[). For further discussion and analysis of flux-based wave decomposition methods, in the context of 
more general approximate Riemann solvers, see pj] ]. While it is difficult to formally prove second order convergence 
for numerical methods since this also requires proving stability, our numerical tests of these methods, and those of 
reference p0| , typically exhibit second order convergence. 



1. Flux-Based Wave Decomposition 

Using Eq. (|54|) to update average grid cell values of the variables q requires knowing flux values at grid cell interfaces. 
The interface flux values are found by solving the following equation, obtained by multiplying Eq. (|54| ) by A(x) on 
the left hand side: 

9 t [F(q)]+A(x)9 x [F(q)] = 0. (55) 

The time derivative of A(x) vanishes because the variables on which A(x) depends (the lapse, the shift, and the 
longitudinal part of the spatial metric) have no fluxes, and are not updated during the transport step. Using Eq. ( p55| ) 
to compute the interface fluxes was originally described by Bona et al. However, it is not clear from || how they 
handled problems in which A(x) varies from cell to cell. 

Eq. flSq ) is a linear advection equation for the flux vector, F(q). As such, flux values at cell interfaces can be updated 
by solving Riemann problems based on decomposing flux differences between adjacent grid cells into eigenvector 
expansions (see H for a discussion of solving Riemann problems for the advection equation) , and including correction 
terms to give second order accuracy. We develop two wave propagation methods based on this idea which we call 
Methods I and II. A wave in this approach is defined as a discontinuity in the flux associated with a certain eigenmodc 
across the characteristic corresponding to that eigenmode. 

We explicitly deal with the fact that the eigenvalues and eigenvectors of the characteristic matrix are varying across 
the grid. The magnitudes of the eigenvalues give the wave speeds and the signs of the eigenvalues give the wave 
directions. In the flux decomposition for Method I, we need to decide if a wave is left-going or right-going at a 
given cell interface. This is determined by the sign of the average of the eigenvalues obtained from the characteristic 
matrices on either side of the interface. If the average eigenvalue for a particular eigenmode is negative, then the 
corresponding eigenvector is evaluated in the cell to the left of the interface. If the average eigenvalue is positive, then 
the eigenvector is evaluated in the cell to the right of the interface. In Method II, the eigenvalues and eigenvectors at 
a cell interface are obtained from the characteristic matrix at the interface, calculated as an average from the adjacent 
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cells. For both methods, waves with zero interface speed still contribute to the flux difference. We can include these 
contributions in either the left-or right-going waves of Method I, as long as we do so consistently. 

In Method I, the flux difference decomposition takes the following form at the interface between cells i and i — 1: 

m M 

F(q i )-F(q i _ 1 ) = Aiq i -A i _ 1 q i _ 1 = ^ 7 f_ 4 rf_ 1 + ^| r f- ( 56 ) 

L=l R.=m+1 

where r are right eigenvectors of the characteristic matrix, and M is the total number of eigenmodes. We denote 
the left-going waves at this interface as W L _j_ = J L _i?f-i, where 1 < L < m. The right-going waves are given by 

W^i = 7^.1^, where m + 1 < R < M. The number of left-going waves, m, can vary from interface to interface 

1 2 1 2 

since the sign of the average eigenvalue can change from cell to cell. The eigenvectors r^_i are evaluated in cell 
Likewise, are evaluated in cell i. The coefficients 7j_ i are obtained by solving Eq. d56|); the subscripts i — | indicate 
interface values. In Method II, the flux difference decomposition at a given interface between cells i and z — 1 is the 
same as Eq. (|56|), except the eigenvectors r f _i of the averaged characteristic matrix A { _i = (Aj_x + Aj)/2 replace 
both rj_i and rj. 

Method I is implemented in the context of the CLAWPACK software package |B2[ . The first order wave propagation 
and second order corrections in both Methods I and II are analogous to Eqs. (18) and (19) of LeVeque's paper on 



standard wave decomposition methods |23|. The updated value of is given by 

(57) 



*=*-^(Ew* 1+ ^)-^( f <+s-f< 



F i± i are flux correction terms which can be reduced near discontinuities by introducing limiter functions. Limiters 
prevent the oscillatory behavior around discontinuities seen with finite difference methods. In the absence of limiting, 
the flux corrections are 



lfc'ii-E" t J-i5EWi' 

\ R L I p=l 



where A^ ±i denote cell-interface speeds. 

Both flux-based wave decomposition methods I and II are successful in giving second order convergent results in 
our numerical calculations. 



2. Finite Difference Method 

To solve Eq. ( |54| ) using a Lax-Wendroff finite difference method, we first perform a second order Taylor expansion 
of q around t: 

q(ar, t + At) = q(.x, t) + At d t q(x, t) + ^At 2 d?q(x, t). (59) 

Observe that 

d t q=-d x [A(x)q\, (60) 

and, taking another time derivative, 

d 2 t q = -d x [A(x)(d t q)} = d x [A(x)d x (A(x)q)}. (61) 

Note that the time derivative of A(x) vanishes as in Eq. (p5|) . Plugging these expressions for dtq and d 2 q into Eq. 
© gives 

q(x,t + At) = q{x,t) - At d x [A(x)q(x,t)} + ^At 2 d x [A(x)d x (A(x)q(x,t))}. (62) 

Making the centered finite difference approximation to the derivatives in Eq. (|62] ) , the updated value of q^ is given by 
At , . , . 1 / At N 2 



q^ = qi - ^(A i+ iq i+ i - A,q j: ) + - I — ) [(Aj + A, ;+ i)(A i+ iq 4+ i - A^) - (A4-1 + Ai)(A,qi - A^iq^i) 



(63) 
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3. Boundary Conditions 



Our numerical methods only require values in one ghost cell at each boundary. We obtain values for all variables in 
the ghost cell by quadratic extrapolation from the three adjacent physical cells. Numerical integration of the constraint 
equations from the last physical cell to the ghost cell by the trapezoidal rule is used to correct the constraint quantities 



D 7 



and (K y y + K z z ) in the ghost cell, with iteration to convergence. 



IX. RESULTS 



A. Initial Conditions 



The initial conditions must satisfy the constraint equations, Eqs. (19) and (|2C|). Since the constraint equations are 
differential equations, they require boundary conditions for their solutions. Different choices of boundary conditions 
correspond to different gauge conditions. We choose symmetric boundary conditions which give flat space between 
two waves. This means that we choose (D y y + D z z ) and (K v v + K z z ) to vanish initially between the waves. 

The variables h xx , K x x , and the combinations (D y y — D z z ) and (K y y — K z z ) are freely specifiable. We normally 
take h xx — 1, K x x — 0, and 



0.25In(£a)=5> 



2 (x - x 0i ) 

7T Wi 



sin [ki(x — xoi) + Si] , 



(64) 



for —Wi < (x — XQi) < Wi, and zero outside that range. The x derivative of Eq. 



1) gives (D y y 



D s *)/2. In 



our standard initial conditions for colliding plane waves, one wave is initially on the left and moving to the right, 
K z z ]i — [{D y v — D z z )/\/h xx ]i. The other wave is initially on the right and moving to the left, with 



with [K y v 



[Ky v 



l(D y v D z 



The parameters for < x < 20 are w^ — 4.0, k% = 1.6, Ai — 0.08, x m = 6.0, 



- K,*] 2 

xq2 = 14.0, S\ = 0, and 82 = n. These initial conditions, depending on the variables, are symmetric (or antisymmetric) 
about x — 10, and symmetry (or antisymmetry) is preserved throughout the evolution. Hence, our figures only show 
the range < x < 10. Since the initial plane waves do not overlap, and (D y y + D z z ) and (K y v + K z z ) vanish at x = 0, 
the initial conditions are two analytic single plane waves of the type described by Misner, Thorne, and Wheeler [p3|. 



0.15 



0.05 




-0.05 



FIG. 1. Initial conditions for derivatives of the transverse metric. The solid line is (D y l 
) / Vhxx- Note that x = 10 is the center of the grid. 



D z 



and the dashed 



line is 2(D y y + D 



Our initial conditions produce large amplitude, nonlinear colliding gravitational plane waves. Our measure of "large 
amplitude" is that h yy and h zz are substantially different from 1 by the time the waves have traversed the grid. It is 
known that nonlinear planewave spacetimes develop a singularity behind the wave |m|,[35|]. For a single plane wave, 
this is only a coordinate singularity, while for colliding plane waves, a physical singularity also develops. The values 
we take for our wave amplitudes are about as large as possible without allowing a singularity to develop during the 
crossing time of the waves. One can get a feel for this value by asking at what amplitude does a singularity develop 
at the left edge of the grid for a single plane wave exiting the right edge? For a single wave as given by Eq. (|3) with 
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the shape specified by our values for iu, and ki, and a flat metric ahead of the wave, the answer is approximately 
0.11. This is an upper limit, however, because the effects of colliding waves add together in a way which is hard to 
estimate. The initial conditions for {D y v ± D z z ) are shown in Fig. [|. 



B. Comparing Evolution Systems 

1. Testing for the Optimal System 

We experiment with several different formulations of the Einstein equations to determine the factors involved in 
improving the global accuracy of ID colliding gr avit ational plane wave calculations. The basic formalisms we test are 



the modified BM, AY, and ADM schemes of Sec. VI. In all of these schemes, using mixed variables rather than lowered 
variables improves accuracy significantly. We also compare alternative ways of handling the redundant variable V x 
in the BM schemes. V x can be left to evolve independently (no-reset BM), or it can be reset periodically to enforce 
the constraint that V x = D y y + D z z (reset BM). Results have been calculated for a range of values of the coefficient 
n of the energy constraint term in the extrinsic curvature evolution equations, from about —0.4 to 1.0, and in some 
cases for values of n > 1. For the ADM and reset BM schemes, the results near and 1 reflect the breakdown of 
hyperbolicity at these values of n. Results are primarily shown for t = 12 since this is the latest time at which the 
physical waves are largely within the grid. 




2 4 6 8 10 



X 

FIG. 2. Evolution of the metric derivatives. The solid line is at t = 8, the dashed line at t = 10, and the dotted line at 
t = 12. (a) (D y y - D,*)/(WQ (b) (D y y + D x *)/(2y/h^0), and (c) D x x /y/h^. 

Fig. ^ shows the evolution of linear combinations of metric derivatives appearing in the eigenmodes from t = 8 
to t = 12, after the physical waves have finished colliding. In these high resolution (4000 cell) calculations, the 
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numerical errors are negligible on the scale of the graph, and we have verified that all the different formulations seem 



to be converging to the same solution. The quantity (D y y — D z 



c ) is shown in Fig. 



The coordinate 



speed of propagation can be read off the graph: it is roughly two units in x for every two units of time until around 
t = 12, when the coordinate speed of light starts to differ significantly from one. Over the same range of times, the 
quantity (K y y — K z z )/2 is within about 3 per cent of — (D y y — D z z )/(2^/h xx ), close to but not identical to what is 
expected from the left-propagating "physical" eigenmode. In Fig. ||b, we see that the steps in (D y y + D z z )/(2y / h xx ) 
are associated with extrema of (D y v - D z z )/(2y/h x ^). At these times, (K y v + K z z )/2 = (D y v + D z z )/(2y/h^) to 
the left of the physical wave. In the vicinity of the physical wave, (K y v + K z z )/2 has step-like features associated 
with steps in {D y v + D z z )/(2y / h xx ), but ascending to the right. In the region between the waves, (K y v + K z z )/2 
is much larger than {D y v + D z z )/ (2y/h xx ) and increases with time. Fig. ||c shows the evolution of D x x / 1 \Jh xx . The 
prominent feature in this figure is a small residual effect (note that the scale of the graph is 1CP 4 ) of the prominent 
feature in K x x shown in Fig. [9] which survives the near cancellation of K x x in the evolution equation of D x x from 
our shift resetting condition (Eq. (|25|)). Since aK x x — /?' = each time the shift is reset, this feature in D x x /\/h xx 
tends to advect along the hypersurface normals. The generation and modification of features in D x x j 1 \Jh xx is due to 
the different evolutions of aK x x and f3' between gauge resettings. The feature in K x x results from our lapse resetting 
condition, Eqs. J22] ) and (23), when a strong imbalance between the transverse D's and K's occurs near the center 
of the grid as the waves collide, creating negative values for S. This in turn causes the limiter to take effect, which 
allows K x x to dip in the negative direction. 

To compare the overall accuracies of different formulations, we present 1-norms of the energy constraint errors in 
Fig. H and 1-norms of errors in D x x /\/h xx in Fig. ^ at t = 12 for 500 cell grids. The constraint errors are predominantly 
errors in the derivatives of the constraint quantities, and are insensitive to errors in the longitudinal variables. For 
each scheme, the 1-norm errors are plotted for a number of values of the energy constraint coefficient, ranging from 
—0.25 to 0.95 at 0.05 increments. Since our ADM scheme is not hyperbolic for n = 0, the transport steps of the ADM 
calculations are solved with the finite difference numerical method, whereas the transport steps of the BM and AY 
calculations use our flux-based wave decomposition methods. The choice of numerical method makes little difference 
to the results. 




-0°25 0.25 0.5 0.75 

energy constraint coefficient 



FIG. 3. 1-norm errors of the energy constraint plotted against the energy constraint coefficient, n, for several different 
formulations of the Einstein equations. Evaluated at t = 12 with a grid resolution of 500 cells. Note that the largest value of 
n plotted is 0.95. Point A is the formulation closest to the standard BM scheme 

Fig. |^ identifies factors which affect accuracy as measured by 1-norm energy constraint errors. It is apparent that 
using mixed variables improves accuracy significantly in the BM formulation. Similar improvements occur in the 
AY and ADM formulations. The 1-norm energy constraint errors for ADM and the ADM- like reset BM schemes 
are almost identical, differing by only 1 to 2 per cent over the range —0.25 < n < 0.85, and are minimized for 
0.25 < n < 0.80. Point "A" on Fig. || is the formulation closest to the BM scheme as implemented in reference [Q. 
The identical formulation using mixed instead of lowered variables decreases the 1-norm energy constraint error 3.1 
times. If the mixed BM system of equations is transformed into an ADM- like scheme by frequently resetting V x , and 
an energy constraint coefficient of 0.5 is used, a 9.3-fold decrease in the 1-norm energy constraint error compared to 
point "A" is obtained. Both the ADM and the reset BM error curves peak at n = 0, and increase rapidly as n — > 1, 
though the increase as n —> 1 for mixed reset BM occurs too close to n — 1 to be apparent in Fig. 0. The rise in 
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energy constraint errors at n = and n = 1 reflects in part the breakdown of hyperbolicity in ADM at these values 
of n. The effects of this breakdown are more severe at n = 1 than at n = because ADM is unstable for n > 1. 
Momentum constraint errors are similar to or smaller than the energy constraint errors. 

The 1-norm energy constraint errors in the no-reset BM schemes vary slowly for all n. These schemes are well- 
behaved for n > 1, and the errors for the mixed version continue to decrease. Despite the breakdown in the AY scheme 
at n = 1, the constraint errors do not increase strongly until n gets close to 1. The AY formulation is well-behaved 
for n > 1. 

Since the true value of D x x /^Jh xx is not known exactly, we must extrapolate to estimate the true value and calculate 
the 1-norm errors shown in Fig. || Assuming quadratic convergence, the error estimate at each grid cell of a 500 cell 
calculation is | times the difference between the 500 and 1000 cell results. For D x x / \/h xx at time t = 12, the 500 
cell errors deviate from quadratic scaling in the region 8 < x < 10, where the feature in D x x /^Jh xx associated with 
the spike in K x x is located. Here, our standard extrapolation underestimates the errors for the ADM and no-reset 
BM formulations, and overestimates the errors for the reset BM formulations, but the effects on Fig. |^ are not very 
significant. 



.x 10 




0.25 0.5 
energy constraint coefficient 

FIG. 4. 1-norm errors of D x x /\/h xx plotted against the energy constraint coefficient, n, for several different formulations of 
the Einstein equations. Evaluated at t = 12 with a grid resolution of 500 cells. Errors are estimated from comparisons with 
1000 cell calculations, assuming quadratic convergence. Note that the largest value of n plotted is 0.95. The legend is the same 
as in Fig. [] except that the AY mixed values are multiplied by 0.045. 

Fig. ^ shows that the same factors which decrease the 1-norm energy constraint errors also decrease the 1-norm 
errors in D x x j \jh xx . Specifically, using mixed variables instead of lowered variables in the no-reset BM formulation 
at n — decreases the 1-norm error in D x x j \Jh xx 1.5 times. The shapes of the curves fall into the same two classes 
as described above for Fig. ^. If one frequently resets V x and sets n = 0.5 in mixed BM, the 1-norm error decreases 
6- fold from point "A" . The errors in D x x / \Jh xx in the ADM and reset BM schemes peak at n = and increase rapidly 
as n — ► 1, again reflecting, in part, the failure of hyperbolicity at these values of n. 

The fact that Fig. || is at all similar to Fig. |3| is because the D x x / y/h xx err ors and the constraint errors behave 
similarly in the region of the physical wave (this is described in detail in Sec. IX B 2). However, there are several 



differences between these two figures. First, the ADM curve in Fig. [| is well above the reset BM curve. Second, the 
reset BM curves have smaller drop-offs from n — in Fig. |[ Third, the decreasing slopes of the no-reset BM curves 
are bigger in Fi g. ^ than in Fig. These differences are due to relatively large formulation-dependent numerical 
errors in D x x j\Jh xx in the region 8 < x < 10, which contribute to the 1-norms. Recall that the evolution of the 
feature in D x x /y/h xx in this region depends mainly on the evolutions of K x x and (3' between gauge resettings. We 
expect the numerical errors for K x x to differ from those for /3', since these two variables evolve by quite different 
equations. Further, we expect these numerical errors to be larger where K x x (and /?') vary rapidly (see Fig. ^|). Small 
differences in these numerical errors among the different formulations result in large differences in numerical errors 
for D x x /y/hxx in this region. For example, the errors in D x x /\fh xx for 8 < x < 10 are larger at n = than at n = 0.5 
for the mixed no-reset BM scheme, explaining the decreasing slopes in Fig. [|, and are larger for ADM than for reset 
BM, explaining the displacement between the mixed ADM and mixed reset BM curves. 

Another difference between Figs. ^ and ^ is that the 1-norm errors of D x x /\/h xx are one to two orders of magnitude 
higher for AY than for the other formulations, whereas the 1-norm energy constraint errors for AY and the other 
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formulations are comparable. For example, the 1-norm error in D x x / 1 \Jh xx at n = 0.5 is about 70 times higher for the 
mixed AY formulation than for the mixed no-reset BM formulation, whereas the 1-norm AY energy constraint error is 
about 2.5 times higher. The subtraction of two numbers with large errors in Eq. ([39]) gives a large error for D x x / \Jh xx 
in the AY formulation. When we introduce substantial variations in h xx in the initial conditions, so that D x x is much 
larger than (D y v + D z z ), and h xx is substantially different from 1, then equally large errors are introduced into all 
the formulations. In our particular calculations, the errors already present in the AY scheme fortuitously cancel the 
introduced errors, resulting in D x x / ' \Jh xx errors for mixed AY and mixed no-reset BM which are within a factor of 2. 

These results demonstrate our ability to significantly increase the accuracy of ID highly nonlinear colliding grav- 
itational plane wave calculations through equation formulation. In particular, the choice of a mixed set of indices 
improves the accuracy of the formulations for all values of the energy constraint coefficient tested. Resetting V x in 
the BM formulations creates ADM-like schemes. This is an advantage for 0.25 < n < 0.80, where ADM is hyperbolic, 
and a disadvantage for n = or n = 1, where ADM is not hyperbolic. 



2. Error Propagation 

In order to understand the variations in accuracy among the formulations, it is instructive to look at how errors 
vary with position, and how they propagate over time. We focus on the energy constraint errors and the errors in 
D x x I \/h xx , since they are representative of errors in the transverse and longitudinal parts of the metric, respectively. 
Momentum constraint errors are comparable to or less than the energy constraint errors. 

Since the principal parts of the energy and momentum constraints are the derivatives of the constraint quantities, 



the constraint errors propagate with the same speeds as the "constraint" eigenmodes (which are given in Sec. VI D 
for BM, AY, and ADM). The constraint error propagation can also be obtained from the evolution equations for the 
energy and momentum constraints, which form their own closed hyperbolic system. The constraint errors propagate 
along the light cones for the no-reset BM and AY formulations. For ADM, the constraint errors propagate at 



vadm = ~P± /7 - — Vl - n. (65) 

We expect the constraint errors to also propagate at vadm for the reset BM scheme. Note that for n < or < n < 1, 
vadm is different from any of the other characteristic speeds of the system. We find that a separation of the constraint 
error speeds from the other characteristic speeds improves accuracy. 

Fig. |]a shows the energy constraint error propagation for reset BM and n — 0. First, notice that the constraint 
errors are large where the physical wave is present, and that at a given location, the error decreases almost to zero 
when the physical wave has passed. However, there is a rapid increase in the energy constraint errors propagating with 
the physical wave. Second, from the graph we see that the waveform of the energy constraint errors propagates at 
roughly unit coordinate speed, which, over the times we are considering, is approximately coordinate light speed. The 
energy constraint errors are predominantly errors in the derivatives of {D y v + D z z ) / '\/h xx . Errors in the propagation 
of (D y v + D z z ) I \Jh xx are largest where its second derivative is largest, at the corners of the steps visible in Fig. 
||b. From the energy constraint equation, the steps in (D y y + D z z )/ \/h xx are associated with large values of the 
physical quantities (D y v — D z z )/\/h xx and (K y v — K z z ). These physical quantities propagate at light speed, and 
constraint errors, once generated, propagate with the velocity of the "constraint" eigenmodes, which is also light speed 
for reset BM with n = 0. Since new errors remain in phase with the propagating old errors, the constraint errors 
are continuously reinforced. Careful comparison of Fig. Hbwith Fig. ||a shows that the constraint error peaks are 
coincident with the corners of the steps in (D y v + D z z )/^/h xx at all three times shown. 

The same argument applies to AY and no-reset BM, since these formulations also have "constraint" mode errors 
propagating at light speed, but Fig. || shows larger errors for ADM and reset BM at n — than for AY and no- 
reset BM. We attribute the larger ADM and reset BM errors to the breakdown of hyperbolicity in ADM at n — 0, 
so that the "constraint" eigenmodes, which propagate at light speed with constant amplitude, interact with the 
"longitudinal" eigenmodes through the constraint quantities. As n — > 0, any errors in the constraint quantities tend 
to produce amplified errors in the longitudinal variables, because n(D x x + Dq) and n(K x x — (3' /a) are the same order 
of magnitude as the constraint quantities in the ADM "longitudinal" eigenmodes. The errors in D x x and K x x then 
feed back to the constraint quantities through the source terms of the evolution equation for (K y y + K z z ). 

Fig. ^b shows a dramatic decrease in both the energy constraint errors and the growth rate of the errors in the 
mixed reset BM formulation when n = 0.5. Furthermore, the full waveform of the energy constraint errors does not 
maintain its shape as it propagates, as does the waveform in Fig. ^|a. By tracking the rightmost bump in this figure, 
one can determine the speed of the errors to be approximately 0.7 times light speed, which agrees with the value 
predicted by Eq. (pa). Because the energy constraint errors lag behind the source of the errors, namely, the steps in 
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(D y v + D z z )/ \/h xx , there is not constant reinforcement and rapid growth of the errors in the region of the physical 
wave. This results in greater overall accuracy for reset BM than for no-reset BM at n = 0.5, as seen in Fig. || 

Fig. [|c shows the energy constraint error propagation for mixed no-reset BM at n — 0.5. In contrast to Fig. |]b, 
the waveform is maintained reasonably well, and the errors grow more rapidly in time. This is because the constraint 
errors travel at light speed; so, as in Fig. ^a, there is a continuous reinforcement of the errors. For no-reset BM at 
n = 0, the curve is practically the same as what we show here at n — 0.5. The error waveform has a smaller growth 
rate than that for reset BM at n = 0, presumably because of the hyperbolicity of the no-reset formulation, as discussed 
earlier. We have also looked at the AY constraint error propagation and confirm that errors propagate at light speed 
for all values of n. Because the constraint errors travel with the physical waves for all n in these formulations, their 
1-norm errors vary slowly with n in Fig. p[ 
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FIG. 5. Energy constraint error propagation using our mixed BM formulation. The solid line is at t = 8, the dashed line 
at t = 10, and the dotted line at t = 12. Evaluated with a grid resolution of 500 cells, for (a) n — 0, with V x resetting, (b) 
n — 0.5, with V x resetting, (c) n = 0.5, with no V x resetting. 

The errors in D x x / ' \Jh xx versus x at t = 12 for the reset BM formulations are shown in Fig. |[ for energy constraint 
coefficients of (a) and (b) 0.5. Fig. ||a shows a localization of the errors in the region of the physical wave (0 < x < 6). 
Fig. ^b shows the dramatic decrease in the errors for < x < 6 when n — 0.5. The mixed no-reset errors in this 
region at n = 0.5 are larger by a factor of about 2 than the mixed reset errors at n = 0.5, and only slightly smaller 
than the mixed no-reset errors at n = 0. The fact that there is a similar reduction in constraint errors and D x x / ' \Jh xx 
errors when going from no-reset BM to reset BM at n = 0.5 suggests that the constraint errors are a major source 
of errors for D x x / \/h xx . From the failure of hyperbolicity in ADM at n = 0, one expects the errors in D x x j \Jh xx to 
increase more rapidly for reset BM than the constraint errors, but we do not see clear evidence for this. The increase 
in errors is about the same going from t = 8 to t = 12, perhaps because the D x x j\fh x ~ x errors have not yet reached 
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their a symptotic limit. The spikey errors at approximately 8.7 < x < 10 in both Figs. ^|a and b are discussed in Sec. 
K B ll 
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FIG. 6. Errors in D x x /\/h xx versus x for our BM mixed and lowered formulations, with V x resetting. Evaluated at t = 12 
with a grid resolution of 500 cells. Errors are estimated from comparisons with 1000 cell calculations, assuming quadratic 
convergence, for (a) n — 0, and (b) n = 0.5. 

The ADM and ADM-like reset BM systems show rapid increases in errors as n — > 1 in Figs. |^ and ^. The 



breakdown of hyperbolicity at n — 1 results in errors in (D y v + D z 



becoming large compared to errors in 



(K y v + K z z ). The contribution to the energy constraint errors from errors in the derivatives of (D y 
increases correspondingly. Since the "constraint" eigenmodes propagate along the hypersurface normals at n = 1, 
the errors in (D y y + D z z )/^/h xx also reinforce errors in variables which propagate along the hypersurface normals, 

namely D x x /\fh^, ft, and Dq. 

Fig. shows the errors in (D y y — D z z ) / (2y/h xx ) as a function of x for three formulations at an n of 0.5. The 
amplitudes and shapes of the errors for the different formulations are roughly the same, because the evolution equations 
for (D y v — D z z ) I '(2y/h xx ) are similar for the different formalisms. The curve shapes do not change significantly when 
one uses n = instead of n — 0.5 because n does not enter into the evolution equations for (D y v — D z z )/ (2\Jh xx ). 
The entire graph shown in Fig. [?] converges quadratically except for the bump around x = 6, which converges linearly. 
This bump is near the trailing edge of the physical wave, where the initial conditions are not smooth enough to give 
second order accuracy. The errors for the formulations which are not shown are similar. 
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FIG. 7. Errors in (D y v — D z z ) / (2\/h xx ) versus x for different formulations. Evaluated at t = 12 with a grid resolution of 
500 cells. Errors are estimated from comparisons with 1000 cell calculations, assuming quadratic convergence, for n = 0.5. 
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Calculations of the distribution and propagation of energy constraint errors and errors in D x x j\Jh xx have illumi- 
nated how resetting V x in the BM formulations to create overall ADM-like evolutions increases or decreases accuracy 
and stability, depending on the value of n. Intermediate values of n separate constraint error speeds from the other 
characteristic speeds of the system, resulting in significant improvements in accuracy. Values of n for which ADM 
hyperbolicity fails result in rapid increases in errors. 

The higher accuracy of reset BM compared to the other formulations when 0.25 < n < 0.80 is seen when h xx — 1 
and D x x — in our initial conditions. When we introduce substantial variations in h xx in the initial conditions, so 
that D x x is at least as large as (D y y + D z z ), and y/h xx varies by a factor of 2 to 3, the accuracy results are dominated 
by the errors directly associated with the variations in h xx . These errors depend primarily on the amplitude of 
D x x / \/h xx , and are largest where the derivative of D x x j \Jh xx is largest. They are similar in all the formulations and 
are independent of n; thus, they tend to equalize the results from the different formulations. 

C. Gauge Conditions 

Our gauge conditions involve both a periodic resetting of the lapse and the shift, and an evolution of the lapse 
and the shift between resettings. The lapse and shift are reset according to Eqs. (g2|) and (p5|), with T = 1, and 
(K x x )t = —0.02. They are reset at constant time intervals rather than after a specific number of time steps, to give 
resolution-independent results. Between resettings, f3, /?', Q, and Dq are advected along hypersurface normals. 

Fig. H shows the behavior of the lapse and the shift at the left boundary of the grid with time. The values of ln(a) 
and (3 are fixed at zero at x = 10, the grid center. ln(a) and (3 vary approximately monotonically from the grid center 
to the edges; therefore, we illustrate how they behave as a function of time by graphing their values at the edge, where 
typically ln(a) and (3 are largest in magnitude. 
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FIG. 8. The shift, /3, and the logarithm of the lapse, tn(a), at the left grid edge plotted against time at intervals At = 0.2. 
The values of f3 and ln(a) are fixed at zero at the grid center. 

In Fig. |8|, the lapse oscillates at the grid edge until about t = 6. From t = 6 to t = 12, ln(a) becomes increasingly 
negative until it starts leveling off around t = 12. The lapse is determined by Eq. (^2|). Since we set the lapse to 
one and its first derivative to zero at the center of the grid, whether the lapse is greater than or less than one at the 
edge is determined by the sign of the second derivative of the lapse with respect to proper distance, the quantity S 
in Eq. (p2|). The large terms in S are typically —K y y K z z and D y y D z z /h xx . These nearly cancel when D y v , D z z , 
K y v , and K z z are dominated by a single propagating "physical" wave pulse. As two such pulses collide, constructive 
interference in D y y and D z z implies destructive interference in K y v and K z z , and vice versa, which causes rather 
large oscillations in the value of the lapse at the edge of the grid. Note that at t = 4, there is maximum constructive 
interference of D y v and D z z . Once the waves have largely passed through each other (t > 6), K y y and K z z are 
large, and D y y and D z z are small in the region between the separating waves, so the second derivative of the lapse 
is strongly negative, and becomes more strongly negative as K y v and K z z become steadily larger. By t = 12, the 
limiter in Eq. ( |22| ) , which is necessary to prevent the lapse from becoming negative at the edge of the grid, largely 
stops further decrease of the lapse at the edge. 

The behavior of the shift in Fig. || reflects the behavior of K x x in Fig. [|, since /3' is proportional to K x x by Eq. 
(p5|). The resetting of the lapse keeps K x x near its target value of —0.02 until t > 10, when the limiter in Eq. ( |2^ ) 
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takes hold, allowing K x x to become steadily more negative. Since the derivative of the shift is proportional to K x x , 
and the shift is clamped to zero at the center of the grid, the shift at the left edge of the grid is positive and starts 
getting steadily larger for t > 10. 
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FIG. 9. K x x versus x at the times indicated in the legend. K x x is driven to a target value of —0.02 by our lapse resetting 
condition. The spike at x ~ 9 develops as the waves collide, due to the limiter in the equation for the lapse. The effects of the 
limiter are felt over a wider range of x at t = 12. 

Recall that a negative K x x , in combination with our shift resetting condition, causes the hypersurface normals to 
point away from the computational grid at the boundaries (see Sec. ^). These resetting gauge conditions are designed 
so that Q and f3 are advected out of the grid, in order to suppress the development of instabilities associated with 
extrapolation at the grid edge. As long as (K x x )t < 0, the shift is positive at the left edge of the grid and the 
advection velocity is negative. However, Fig. [l0| shows that when (K x x )t = +0.02, so that the shift is negative at the 
left edge, the solution becomes unstable at the grid edge once the physical wave reaches the edge. 
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FIG. 10. Development of an instability at the left boundary in the evolution of D x x / \fh xx when the K x x target value is 
+0.02 in our lapse resetting condition. The instability is substantially stronger with the higher resolution grid. Note that the 
physical wave reaches the left edge at t = 10. 



D. Boundary Conditions 



Fig. ^ compares two methods for implementing boundary conditions: quadratic extrapolation of all the variables 



and (K y v + K z z ) to insure that the energy 



with and without correcting the extrapolated values of (D y v + D z z )/y 
and momentum constraint equations are satisfied at the boundaries. Failure to strictly enforce the constraints at the 
edges of the grid results in significant constraint errors propagating over much of the grid by t = 8. 
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FIG. 11. Energy constraint error for our standard boundary condition, which corrects quadratic extrapolation of the variables 
at the boundaries using the constraint equations (solid line), versus quadratic extrapolation only (dashed line). Evaluated at 
t = 8 with a grid resolution of 500 cells using our mixed BM formulation with Vx resetting and n = 0.5. 

The jitter at the left edge of Fig. ^], which graphs the errors in D x x / ' \Jh xx versus x, is due to numerical errors 
and the use of quadratic extrapolation in our boundary conditions. It is most prominent when the physical wave is 
crossing the boundary, but it is not unstable as long as the hypersurface normal at the boundary does not point into 
the grid (see Sec. |v|). We will explore ways of reducing this jitter either through the use of limiters near the boundary 
and/or improved extrapolation methods. 

E. Numerical Methods 



TABLE I. Convergence results for errors in the metric derivative variables. Evaluated at t = 12 using our mixed BM 
formalism, resetting V x , with n — 0.5, and flux-based wave decomposition method I. 
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(500 - 1000)/(1000 - 2000) a 
(1000 - 2000)/(2000 - 4000) b 




3.98 




3.99 


5.09 




3.97 




3.99 


4.38 



a The 1-norm of the difference between the 500 and 1000 cell calculations is divided by the 1-norm of the difference between 
the 1000 and 2000 cell calculations. 

b The 1-norm of the difference between the 1000 and 2000 cell calculations is divided by the 1-norm of the difference between 
the 2000 and 4000 cell calculations. 



TABLE II. Convergence results for the energy and momentum constraint errors. Evaluated at t — 12 using our mixed BM 
formalism, resetting V x , with n = 0.5, and flux-based wave decomposition method I. 
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500/1000 a 


4.02 


3.97 


1000/2000 b 


4.00 


3.98 


2000/4000 c 


3.99 


3.99 



a The 1-norm of the 500 cell calculation is divided by the 1-norm of the 1000 cell calculation. 
b The 1-norm of the 1000 cell calculation is divided by the 1-norm of the 2000 cell calculation. 
c The 1-norm of the 2000 cell calculation is divided by the 1-norm of the 4000 cell calculation. 
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We compare convergence results from our two flux-based wave decomposition schemes with those from the tra- 
ditional finite difference approach. Tables [| and [n] show quadratic convergence for errors in the metric derivative 
variables and in the momentum and energy constraints at t — 1 2 using our optimal evolution scheme and flux-based 
wave decomposition method I. Quadratic convergence for the metric derivative variables is tested by comparing I- 
norms of differences between results evaluated at 500 and fOOO cell resolutions, fOOO and 2000 cell resolutions, and 
2000 and 4000 cell resolutions. Convergence is calculated in this way because true values are unavailable for the metric 
derivative variables. Constraint errors can be calculated directly so quadratic convergence is determined simply by 
taking the ratios of results from grid resolutions that differ by a factor of two. There is very little difference in the 
results of the three numerical methods tested. Further, the f-norm results from the three numerical methods all 
converge to the same answer. 

Tables | and || show that our flux-based wave decomposition methods are second order convergent when the 
characteristic matrix A(x) has significant spatial variation, because of the spatial dependence of the lapse and the 
shift (sec Eq. ( p7[ ) and Fig. ||). However, in these calculations the eigenvectors of A(x), which depend only on 
\/h xx , are roughly constant. We have tested convergence of both our flux-based wave decomposition methods when 
the eigenvectors of A(x) vary rapidly, by introducing a 3-fold variation in \Jh xx . Our results are still second order 
convergent. 



X. DISCUSSION 



A. Summary 

We have identified ways of improving the accuracy and stability of f D nonlinear colliding gravitational plane wave 
calculations through an in-depth study of equation formulations, dynamic gauge conditions, boundary conditions, and 
numerical methods. Three issues stand out in our study of equation formulations. The first issue which improves the 
accuracy of all the formulations tested is raising an index in the metric derivative and extrinsic curvature variables. 
The separation of constraint and physical behavior in our calculations is much simpler with mixed variables than 
with lowered variables. Since the dominant feature of plane wave solutions is the physical wave, and since the mixed 
variables have a simple linear relation to the "physical" eigenmodes, the mixed form gives an advantage in accuracy. 
In addition, the variables in the source terms of the evolution equations are in mixed form, so a cleaner system of 
equations results. 

Second, it is advantageous if the "constraint" eigenmodes, which are constraint violating, propagate at different 
speeds from other features in the solution, so that nothing in the actual solution is constantly in step with the constraint 
errors. Our gauge conditions result in features which propagate on the hypersurface normals as well as on the light 
cones. With a multiple between 0.25 and 0.80 of the energy constraint equation added to the evolution equations for 
the extrinsic curvature in the ADM and reset BM schemes, the "constraint" eigenmodes propagate neither on the light 
cones, nor on the hypersurface normals, and the growth rates of constraint errors and errors in D x x / ' \Jh xx significantly 
decrease. However, when the amplitude of D x x is large compared to {D y y + D z z ), formulation-independent errors 
associated with the derivatives of D x x j\Jh xx dominate the overall errors for 0.25 < n < 0.80. 

Third, we find that when hyperbolicity fails in ADM, for values of the energy constraint coefficient equal to or 1, 
the constraint errors and errors in D x x j \J h xx increase much more quickly in both ADM and reset BM than they do 
in no-reset BM. At n = 0, the "constraint" and "longitudinal" eigenmodes are no longer independent, causing a rapid 
growth of errors which travel with the physical waves. These errors decrease and stabilize, however, when the waves 
exit the grid. At n = 1, the two "constraint" eigenmodes are not independent, causing errors in {D y v + D z z ) / ' \Jh xx 
to increase rapidly. This results in large errors in the energy constraint and in variables which propagate along the 
hypersurface normals. 

The key to our approach to dynamic gauge conditions for hyperbolic calculations is to maintain a simple diago- 
nalizable hyperbolic evolution during each time step, while allowing periodic, flexible resetting of the lapse and the 
shift between time steps. By resetting the lapse and the shift, we can control the coordinate system in a dynamic 
and (in principle) arbitrary way, unconstrained by the need to maintain hyperbolicity. Our gauge resetting conditions 
control the longitudinal components of the extrinsic curvature and the spatial metric so as to prevent pathologies and 
strong gradients from developing in the hypersurfaces and spatial coordinates. Further, our gauge resetting conditions 
cause the hypersurface normals to point away from the grid at the edges. This helps to suppress the development of 
instabilities at the boundaries which are associated with features advecting along the hypersurface normals. 

A careful study of boundary conditions has led us to the conclusion that it is incorrect to impose outgoing wave 
boundary conditions based on the eigenmodes of the hyperbolic decomposition. A substantial contribution from 
the incoming "constraint" eigenmodes is necessary to satisfy the constraint equations at the boundaries. In the 
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presence of nonlinearities, the interaction of the "physical" eigenmodes with the incoming "constraint" eigenmodes 
generates an admixture of incoming and outgoing "physical" eigenmodes. The right and left-going gravitational wave 
amplitudes can be determined by projecting the Weyl tensor onto a null tetrad, and our numerical results indicate that 
the gravitational waves after the collision are purely outgoing according to this definition. Our boundary condition 
procedure consists of an accurate calculation of the incoming eigenmodes of the characteristic matrix at the boundaries 
by quadratic extrapolation of all the variables, and correction of the ghost cell values using the energy and momentum 
constraint equations. This procedure, in combination with our lapse and shift resetting conditions which insure that 
the hypersurface normals at the boundaries do not point into the grid, has given stable, accurate, and second order 
convergent results. 

Finally, we have developed flux splitting numerical methods for solving hyperbolic formulations of the Einstein 
equations which are second order accurate for smooth solutions, even when the eigenvalues and eigenvectors of the 
characteristic matrix are spatially varying. These methods are based on decomposing flux differences between adjacent 
grid cells into linear combinations of the eigenvectors of the characteristic matrix. We show that these methods are 
formally second order accurate and, in practice, second order convergent. 

B. Relevance of Results and Future Directions 

Our results suggest that hyperbolicity can be a useful guide to picking equation formulations for numerical integra- 
tion of the Einstein equations. When eigenvectors are incomplete, or some eigenvectors are nearly linearly dependent, 
some solutions to the equations will tend to grow without bound or by large factors. In 2D and 3D, the ADM and 
ADM-like equation formulations for non-diagonal metrics fail to have a complete set of characteristic eigenvectors for 
any value of the energy constraint coefficient. This may be at the root of some of the instabilities seen in 3D ADM 
codes. 

However, hyperbolicity, which involves only the principal terms in the equations, is not the whole story. The 
constraint equations relate derivatives of some quantities (the constraint quantities) to nonlinear terms involving 
additional quantities, so the actual evolution of the constraint quantities for constraint-satisfying solutions may be 
very different from the evolution implied by the "constraint" eigenmodes, which satisfy equations without source 
terms. It seems to be advantageous if the speeds of the short-wavelength errors in the constraint quantities, which 
are the "constraint" mode eigenvalues, differ substantially from the propagation speeds of the major features in the 
constraint quantities. In ID, both the ADM and the no-reset BM formulations are safely hyperbolic for n ss 0.5, but 
the former is substantially more accurate because no-reset BM has the same speeds for the "constraint" eigenmodes 
and the features in the constraint quantities, whereas ADM has different speeds. How much of an advantage it is to 
have different speeds depends on the dominant source of errors. In a gauge where h xx is close to 1, the dominant 
numerical errors in the constraint quantities are associated with features in the same quantities, but in a gauge where 
the dominant numerical errors are associated with features in h xx , it does not make much difference whether the 
"constraint" eigenmodes and features in the constraint quantities have the same velocities. In generic black hole 
spacetimes, it is probably not possible to find a gauge where the physical waves dominate the metric; therefore, 
the separation of speeds may not make much difference in accuracy for 2D and 3D calculations of greatest physical 
interest. 

The eigenmode decomposition in higher dimensions depends on a chosen direction of propagation, as do the com- 
binations of the actual variables which can be identified as "longitudinal", "constraint", and "physical" quantities. 
In addition, putting the equations into first-order form requires a choice of derivative ordering, which has important 
effects on the hyperbolicity of the system. A "directional splitting" approach to solving the transport steps, in which 
d t <\+ dkF k — is solved separately for each coordinate direction k, does, at least for some choices of derivative order- 
ing, allow the identification of "longitudinal" , "constraint" , and "physical" eigenmodes for one coordinate direction 
at a time. The decomposition involves projecting the and K^s perpendicular to and into the constant-x fe 

surfaces. The eigenvectors can be constructed explicitly. 

Using the mixed coordinate components Dki and Ki° as variables when the metric is not diagonal gives complicated 
source terms in the D^i 3 evolution equations, because the D^i 3 are no longer pure derivatives. Another possibility is 
to take as the variables the components of an orthonormal triad and the projection of the extrinsic curvature onto 
the triad (or the Ashtekar variables Jl5| - |l7| ). The triad formalism, because the extrinsic curvature tensor projected 
on the triad is symmetric, has fewer variables than the mixed coordinate component formalism. In neither case would 
the variables be simply related to the eigenvectors of the characteristic matrix when the metric is non-diagonal or the 
triad vectors are not along the coordinate directions; thus, there is no obvious advantage to trying to generalize our 
mixed variables in ID to generic 2D and 3D calculations. However, we plan to explore these questions further. 

Our approach to dynamic gauge conditions, namely, implementing a simple hyperbolic evolution during each time 
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step, then resetting the lapse and the shift between time steps, is general. A hyperbolic evolution during each time 
step allows the use of flux-based wave decomposition numerical methods. A long-term strictly hyperbolic evolution, 
however, is destroyed in all our hyperbolic formulations when we reset the lapse and shift. Despite this, resetting the 
gauge can be used to increase the accuracy and stability of our ID calculations. Achieving the same goal in 2D or 
3D calculations will be more difficult. There is no one longitudinal direction, and there is not enough gauge freedom 
to control the longitudinal components of the metric and extrinsic curvature for all directions. What gauge resetting 
conditions are most effective in 2D and 3D remains an open question. 

Our boundary condition results indicate that simple outgoing wave boundary conditions based on the eigenmodes 
of the characteristic matrix are not valid in general and are inconsistent with the constraint equations. On the other 
hand, extrapolation (in particular, quadratic extrapolation) is a numerically dangerous procedure. Although one can 
use constraint extrapolation to constrain certain linear combinations of the variables, this technique is likely to be less 
effective in higher dimensions than in ID because the number of variables increases more rapidly than the number of 
constraints. Using an outgoing wave boundary condition based on the Weyl tensor to further constrain the variables 
is problematic in the general case. While the peeling theorem |3(| in an asymptotically flat spacetime does show that 
the incoming wave projection of the Weyl tensor falls off much more rapidly than the outgoing wave projection in 
the wave zone, one cannot assume that the outgoing wave dominates in 3D numerical relativity calculations, which 
typically have to be truncated at best in the inner part of the wave zone. The ingoing part of the Weyl tensor contains 
non-radiative quasi-static quadrupole moment contributions which cannot be assumed to vanish. 

The flux-based wave decomposition numerical methods we have presented for hyperbolic formulations of the Einstein 
equations are completely generalizable, and should prove useful for calculations in black hole and Brill wave spacetimcs 
where the eigenvectors and eigenvalues have significant spatial variation. We hope to extend our exploration of 
hyperbolic methods to include the use of limiters and upwind differencing. Limiters will be used to suppress short 
wavelength numerical instabilities. With upwind differencing, one does not need ghost cells at the apparent horizon 
boundary of an excised black hole, where the eigenmodes are purely ingoing. 

In conclusion, we have developed basic methodologies for hyperbolic formulations of the Einstein equations, which 
improve accuracy and stability in ID, and which we think merit further exploration in the context of numerical 
relativity calculations of more substantial physical interest. 
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APPENDIX A: SECOND ORDER ACCURACY OF FLUX-BASED WAVE DECOMPOSITION 



In the wave propagation stage of solving for the time evolution of a hyperbolic system, the equations being solved 
are exact conservation laws (Eq. (Q)). The integral form of the conservation law can be used to update the average 
q's in the ith cell, 



1 

A~x 



_n+l 



1 



dx = q? - -r- 

Ax 



F„_ i dt 



(Al) 



where i ± | denotes values on the cell boundaries. We base our discussions of accuracy on Taylor series expansions in 
both t and x, assuming that the cell size Ax and the time step At = t n+ \ — t n are the same order. Time and spatial 
derivatives are related by the wave propagation equation, Eq. (^4|), and F is assumed to depend on q and x. Second 
order accuracy means that the error in one time step in q™ +1 decreases faster than (At) 2 as At and Ax go to zero. 
Expanding F in a Taylor series in time, 
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In order that q™ be accurate to second order while neglecting the last term in Eq. (A2), dfF must be continuous 



in x. From the evolution Eq. (|55| ) for F, assuming the dependence of F on q is analytic, second order continuity in 
dfF is guaranteed by second order continuity in d 2 F. The resulting equation is 



t n+ l 



F J4 - 1 dt rs A* F r 



-At 2 A i± i(d x FY 



(A3) 



Let Fi be the cell-centered values of F, which differ from F(qj, Xi) by a term of order Ax 2 d 2 Fi. Then expanding 
F in a Taylor series in x about the cell-center value, 



(^F)«. 



(A4) 



Since the second-order error terms are the same for F i+ i and F^i, and change the same way when Fj and Fj+i 
are replaced by F(q i; Xi) and F(q i+1 , x i+1 ), they cancel in the difference of the flux integrals in Eq. (Al) and can be 
omitted, if d 2 F is continuous. 



To obtain a second-order accurate contribution to q" from the second term in Eq. (A3), it is sufficient that A i± i 



and (d x F)™ ± A be accurate to first order. In Method II, A i± i is approximated to first order by ^(A^ + A i±1 ). A 
Taylor expansion of (d x F) n +L gives, to first order, 
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(A5) 



Flux-based wave decomposition consists of decomposing the flux differences AF" t into a sum over eigenvectors 

*+ 2 

of the characteristic matrix, and similarly for AF"_i. Different versions evaluate these eigenvectors at different 

1 2 

locations. From the point of view of smooth solutions, the preferred location is to base the decomposition of A i± i 
at the respective cell boundaries. The decomposition of AF" +i into a sum of right eigenvectors of the matrix A. i+ i 
can be written 



AF" , = K i+i T i+ i, 



(A6) 



where R,, i is a matrix whose columns are the right eigenvectors r p , , , and the column vector F, • , i consists of the 

l +2 ° ° l +2 5 

coefficients of the decomposition, 7 P ,i ■ Then 



A i+| AF« ,=R i+i A i+ ir 



(A7) 



where A i+ i is the diagonal matrix of eigenvalues of A i+ i. Let A? i be the pth eigenvalue, or wavespeed, and let 



W p ! = 7 ? ! r p ! , where r p ' x is the pth eigenvector, ie., the pth column of R i+ i . Equation (A7) can be written 

A ^ AF iu = E A ^ w ^- ( A8 ) 



Combining Eqs. (|A]|-([Aj) and (|A|), and noting that AF^ +i = J2 P w f+i , gives 
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(A9) 



and a similar expression for the flux integral at i — i. In Eq. ( [Al l ), t hese give for q™ +1 an expression which is easily 
shown t o be equivalent to the substitution of Eqs. ( fig ) into Eq. (|57|). Thus, we have shown that Method II of Sec. 
VIII B 1 is second-order accurate for d 2 F continuous. 



Method I is based on expressions similar to Eq. (A9) for the flux integrals, but the eigenvectors in W i+ i, for 



instance, are derived from Aj for the left-going eig enmodes at the i + ^ interface and derived from Aj+i for the 
right-going eigenmodes. In the last term in Eq. (|A9|), the eigenvalues X p ' x are approximated to first-order accuracy 

I 1 1+ 2 

as the average of the eigenvalues associated with the adjacent cells, but the W j+ i are only zeroth-order accurate. 
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The first-order corrections to the eigenvectors are the same at the i ± i interfaces, if the same eigenmodes are left and 
right-going at each interface and the first derivatives of the eigenvectors are continuous. (Recall that in Method I, 
the sign of the averaged eigenvalue at a given interface determines the direction of the eigenmode at that interface.) 
The second condition normally follows from continuity of d x A., but may fail if some of the eigenmodes are nearly 
degenerate. With this qualification, the first-order corrections in this term will cancel between the two interfaces, as 
long as the same eigenmodes are left- and right-going at each interface. When the averaged eigenvalues at the cell 
interfaces do change sign, continuity of the first derivatives of th e eigenvalues implied by continuity of d x A. means 
X p i+L and A*_ i are both of order Ax. The last term in Eq. ( |A9|) , which becomes the last term in Eq. (f38|), is then 

second-order accurate by itself. The terms first-order in At in Eqs. ( |57| ) and ( |58| ) will differ in second order from their 
values in Method II, but combine by construction to give 
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and 
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(All) 



The result for q™ +1 is the same as Method II through second order, as long as the eigenvectors are smooth as discussed 
above. While eigenmodes are degenerate for many hyperbolic formulations of the Einstein equations, the eigenvectors 
can be chosen with the required smoothness. 

Both methods, when the smoothness conditions are satisfied, are equivalent through second order to the Lax- 



Wendroff finite difference scheme presented in Sec. VIII B 2 
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